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We numerically investigate the motion of active artificial microswimmers diffusing in a fuel con¬ 
centration gradient. We observe that, in the steady state, their probability density accumulates in 
the low-concentration regions, whereas a tagged swimmer drifts with velocity depending in modulus 
and orientation on how the concentration gradient affects the self-propulsion mechanism. Under 
most experimentally accessible conditions, the particle drifts toward the high-concentration regions 
(pseudo-chemotactic drift). A correct interpretation of experimental data must account for such an 
“anti-Fickian” behavior. 

PACS numbers: 


I. INTRODUCTION 

Chemotaxis, defined as the movement of motile cells or 
organisms in response to a chemical gradient, is a well- 
studied phenomenon [l[. Bacteria and other single- or 
multi-cellular organisms propel themselves up or down 
the concentration gradient of a particular substance in 
their search for nutrients or to avoid antagonists. In¬ 
spired by chemotaxis in biology, researchers synthesized 
artificial microswimmers that can move in response 
to a chemical stimulus [J, Q- They showed that Janus 
particles (JP), in the form of two-faced Au-Pt colloidal 
rods that catalyze hydrogen peroxide redox, are attracted 
by a hydrogen peroxide source. Under such conditions, 
JP’s act as molecular “robots” and can thus be employed 
in practical ^plications, such as the design of new intel¬ 
ligent drugs Q . More sophisticated chemical robots have 
been proposed that utilize artificial chemotaxis to navi¬ 
gate autonomously Q- 

The simplest and, possibly, best established model of 
self-propulsion is encoded by the Langevin equations 

na 

X = noCOS(/)-|- a/Ad V = Vq Slli (j)-\- ^/I^ {t), 

^ ( 1 ) 

where r = (x, y) are the coordinates of the swimmer in 
the plane, uq its self-propulsion speed, and an orienta¬ 
tional diffusion constant, whose reciprocal, r^, quantifies 
the time-persistency of the particle’s Brownian motion. 
Here, ^i(t), with i = x,y,(j), are zero-mean and delta- 
correlated Gaussian noises with (Ci(i)^j(O)) = 2SijS{t). 
For long observation times t, with t T,f,, or lengths 

I, with I = Vqt^, the effective diffusion constant, 

J, defined by the asymptotic law limi_,.oo(r^(t)) = 

[ 13 , is U = Dq + Ds, where Dq is due to thermal fluctua¬ 
tions in the suspension, and Dg = Vq/2D^ is a (typically) 
much larger self-propulsion term, which depends on the 
chemical composition of the suspension itself. 




FIG. 1: (Color online) Chemical reactor with a stationary 
fuel concentration gradient, p{x) (see text). A Janus particle 
injected in the middle, (a), tends to drift to the right even 
if its probability density, Po{x), peaks on the left (b). The 
data in (b) are for ao(a;) oc p{x) with Do — 0.01, yv = 1, <5^ = 
1, = 0 and Dij, = 0.1 


Let us consider now a chemical reactor consisting of a 
narrow, straight channel of length L oriented along the 
x-axis, and a free JP moving in it (Fig. [T|) Q. A con¬ 
stant concentration gradient of the chemical that fuels 
the particle’s self-propulsion is maintained by connecting 
the channel to two reservoirs in thermal equilibrium with 
concentrations pL < Pr- The chemical concentration in 
the channel, p{x), will then grow linearly from left to 
right. At the channel ends, x = ±L/2, two porous mem¬ 
branes allow the chemical flow in and out, but prevent the 
JP from escaping into the reservoirs. We speculate, based 
on experimental observations that both vq and 

(and therefore Dg) may depend on p{x) to some un¬ 
specified extent. We only assume that both vo{x) and 
D^{x) are non-decreasing functions of the channel coor- 
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FIG. 2: (Color online) Janus particle in a concentration gra¬ 
dient, Eqs. with rj^ = I, = 2 and different and 
(a) Po{x) for Do = 0.01. Data for Do = 0.01 and different 
are plotted in the inset; (b) Nr/Nl and Tl/Tr (inset) vs 
Do; and (c) e vs Do, Eq. (0. Other simulation parameters: 
Vo = 1, Dcf, = 1, L = lOOZ^ and channel width yi = 5. The 
solid curves are the analytical predictions based on Eqs. 

with a = 1/2 and a — 1, respectively, for <5,^ = 0 and 
= 0 . 


dinate x. Indeed, a higher fuel concentration typically 
enhances active Brownian motion. For this reason, the 
right and left channel endpoints are termed, respectively, 
hot and cold. We then ask ourselves two closely related 
questions. Upon injecting the JP at the center of the 
channel, a; = 0: (1) Which containment membrane is the 
JP more likely to hit first? (2) On which side of the 
channel is it going to sojourn the most time? 

This might sound paradoxical, but we came to the 
conclusion that the injected JP is finally attracted to¬ 
ward the left (cold) exit, even if, immediately after injec¬ 
tion, it may drift to either direction, depending on the 
x-dependence of the propulsion parameters vg and r^. 
For the most common case when the cc-dependence of 
is much weaker than Dg [1, [l3 - [l^ , the injected particle 
points decidedly to the right (hot) exit. Reconciling these 
seemingly conflicting mechanisms is of paramount impor¬ 
tance to control the chemotaxis of artificial microswim¬ 
mers as opposed to bacterial chemotaxis [13, El]. To 
avoid misunderstandings we remark that the adjectives 
hot and cold refer here to the regions in the reactor where 
the effective swimmer diffusion due to the selfpropulsion, 
Ds(x), is the highest and lowest, respectively. The ther¬ 
mal diffusion, Dg, is assumed to be x-independent, which 
means that thermal gradients do not enter our analysis. 
Accordingly, in the absence of thermal gradients and for 
low fuel concentrations, additional transport contribu¬ 
tions due to hydrodynamic effects in the suspension fluid 
can be safely neglected. 

This paper is organized as follows. In Sec. |TT1 we 
present numerical results for the “splitting probabilities” 
that a JP, injected at the center of the channel, exits it 
through the right or left end and the corresponding mean 
first-exit times. The particle clearly undergoes a tran¬ 
sient drift toward the hot end of the channel, whereas its 
stationary distribution tends to accumulate at the oppo¬ 
site end. In Sec. uni we interpret our data by means of a 
phenomenological ID Langevin equation that describes 
the diffusion of a Brownian particle with the spatial de¬ 
pendent diffusion coefficient Dgix). The spatial depen¬ 
dence of Ds{x) generates the drift term here detected as 
a transient drift. Finally in Sec. IIVI we discuss the im¬ 
plications of our findings in the interpretation of recent 
experiments on the diffusion of JP’s in concentration gra¬ 
dients. 


II. NUMERICAL RESULTS 

Our answers to questions (1) and (2) are based on the 
simulation data reported in Fig. [2l As a study case, we 
considered the x-dependent self-propulsion parameters, 

Vg{x) = uo(l + 5yx/L)^'‘, D^{x) = + 5^x/L)^*,{2) 

where Sy = Avg/vg and = AD^/D^ are both non¬ 
negative, and from now on, vg and are shorthands 
for uo(0) and D^/O) at the injection point. We also set 
riy = 1 and = 2, so that for Sy = S^ the self-propulsion 
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diffusion term, Db{x) = Vq{x)/2D^{x), is a;-independent, 
i.e., Db{x) = Dg. That is, the JP is expected to dif¬ 
fuse uniformly along the channel. On the contrary, we 
observed that the stationary probability density function 
(p.d.f.), Po{x), of a single JP in such a close-ended chan¬ 
nel tends to accumulate against the left exit, as displayed 
in Fig. mja). This effect is the strongest as Sy is increased 
at 6^ = 0. Vice versa, as is raised and Sy lowered, 
Po{x) tends to flatten out. For 6y = 0, no matter what 
6^, Po{x) = L~^. Actually, the x-dependence of D^{x) 
seems not to sensibly affect ^ 0 ( 2 ;) for any choice of 
[Fig. HJa), inset]. In conclusion, to answer question (2), 
the injected JP tends to dwell where vo{x) is the lowest, 
as suggested in Ref. that is by the cold extremity 
of the channel (reverse chemotaxis). 

When one looks at the transient dynamics immedi¬ 
ately following the particle injection, a surprising out¬ 
come appears. We injected the particle at x = 0 and 
clocked the time it takes to hit either the right or left 
containment membrane. We repeated this numerical ex¬ 
periment N = 10® times and determined the probability 
the particle first reached the right or left exit, 
and the corresponding mean-first-passage times (MFPT), 
from 0 to ±Ll2. The ratios Nr/Nr and Tr/Tr 
are plotted, respectively, in the main panel and the in¬ 
set of Fig. Iljb). In the regime of low thermal noise, 
Dq Ds, we obtained distinct results, depending on 
which x-dependence is stronger, vo(x) or D^{x). [Note 
that we used the same rjy and as in the main panel (a) 
for Pq (x) .] In the first case, the particle tries to leave the 
channel through the right exit and, accordingly, Tr > Tr. 
Vice versa, on suppressing the x-dependence of vo(x), 
while leaving D^{x) unchanged, the particle directs itself 
preferably toward to left exit and Tr < Tr. 

This means that for Sy ^ the injected particle ini¬ 
tially drifts up the p(x) gradient (pseudo-chemotaxis), 
at odds with Ref. [l3|. Only when the increase of xo(x) 
along the channel is accompanied by a suitably stronger 
increase of its orientational rate, the injected parti¬ 
cle drifts immediately down the gradient, in agreement 
with Tick’s law for ordinary Brownian motion. Magni¬ 
tude and orientation of the transient drift are character¬ 
ized in the forthcoming section by means of the unique 
rectification factor e. This result is remarkable because 
Po(x) tends to accumulate in any case around the con¬ 
centration minima. This behavior is reminiscent of the 
“drift without current” effect experimentally observed in 
Ref. [l^,[13 and numerically investigated in Refs. [2l| for 
thermal Brownian motion in confined geometries. How¬ 
ever, the magnitude of the phenomenon reported here is 
significantly larger and more easily accessible to experi¬ 
mental demonstration. 


III. PHENOMENOLOGICAL ANALYSIS 

An analytical treatment of the model of Eqs. (|T|) is vi¬ 
able in two limiting cases, i.e., 5^ = 0, Sy > 0, and Sy = 0, 





FIG. 3: (Color online) Channel diffusion for r/y = 1 and S^ = 
0: (a) Po{x) for = 0.1, ao = 1.5, and different Sv and Do 
(see legends). Inset: Po{x) for Do = 0.02, D^ = 0.1, = 1 

and different rjy-, (b) Nr/Nr vs. Do for =2/3 and different 
and (c) Tr/Tr vs Do for Sy = 2/3 and different D^. 
Other simulation parameters are: vq = 1.5, L = lOOlrf, and 
channel width yi = 5. Inset in panel (b): Nr/Nr vs D^j, for 
Do = 0.03 and L = 100 (squares) and lOOZ,^ (circles). Inset 
in panel (c): Tr/Tr vs L/l^ for Do = 0.03, and different D^. 
The remaining simulation parameters are as in the relevant 
main panel. Dashed and solid curves are the corresponding 
analytical predictions based on Eqs. (O-® for Q = 1/2. 
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(5^ > 0. For this purpose we implemented the approach 
of Ref. [23 to reduce the fully three-dimensional dynam¬ 
ics of Eq. ([TJ to the more tractable ID phenomenological 
diffusion law, 

X = aD'^(x) + \/Da(x) ^(t), (3) 


fij S,f, = 0, Sy > 0. In view of the third equation ([T]), 
we know that cos (p(t) behaves like a (non-Gaussian) col¬ 
ored noise with an asymptotic autocorrelation function 
(cos0(t)cos(/)(O)) ~ for t ^ [ 2 ^. The JP 

diffusion coefficient at x can thus be derived from Kubo’s 
formula [H, 


where the prime sign denotes an x derivation and (i) 
a = 1/2 and Di/ 2 (x) = Do + Vq(x)/2D^, for = 0, and 
(ii) a = 1 and Di{x) = Do+Vq/2D^{x), for 5y = 0. Here, 
the multiplicative noise term has to be handled according 
to Ito’s prescription [T^ and ^(t) is defined like the noises 
of Eq. dl]). Note that the Eq. ([3]) can be rewritten as 
X = {x) o ^{t), with o denoting the Stratonovitch 

or anti-Ito prescription, respectively, in case (i) and (ii). 
The corresponding Eokker-Planck equation (PPE) is 


-P(x t) = — 
dt ^ ^ > dx 


-Va{x) -b ^Da{x) 
OX 


PM ( 4 ) 


d r ^ 


with Va{x) = aD'{x) for the appropriate value of a [13, 
[ 2 ^. The stationary p.d.f. for zero net current, jo = 
j(x,t) = 0, reads 


Po{x) = livn P[x,t) = N/Da{xY (5) 


D = Do-b lim / v^{x) {cos 4>{t) cos 4){0))dt = Di/ 2 {x), 
as anticipated in Eq. @. 

The drift velocity, Va{x), of a JP with an x- 
independent self-propulsion time constant, amounts 
to the average of vo{x) over its persistence length l^{x) = 
vo{x)t^, i.e.. 


- 2^o(^)^o(a;)T0, 

( 8 ) 

hence Va = aD'^{x) as in Eq. dH) with a = 1/2. 

(ii) Srf, > 0, Sy = 0. Calculating D{x) in this case is 
straightforward. The EPE corresponding to the first and 
third Langevin equations o, 


«a(x) = - 




-VQiX- — 


^p= 

dt 


d 

-v„co.4>-g^ + D„^ + D4r)^ 


where Af is a normalization constant. In particular, for 
0 = 1, i.e., a;-independent vo, Po{x) = 1/L. The exten¬ 
sion to cases with jo ^ 0 is straightforward. 

Regarding the transient statistics of a particle injected 
at the center of the channel, x = 0, a simple “splitting 
probability” calculation leads to 

Nr _ llL/2lP’MPoix)]~^dx 
^ I-i%[DMPo{x)]-^dx 

with Nr + Nr = N. Analogously, for the MEPT’s in a 
channel with absorbing endpoints, we have [l3j | 

TL{Sa) = {T{x)){0,L/2) - {T{x)){-L/ 2,L/2) = TR{-Sa), 

( 7 ) 

where i5i/2 = Sy and =6^, 

/ X pL J'2 

dztjja{z)/Da{z) / dy/ijaiy), 

■L12 J z 

and (...)(a,b) = Ja{---)dx/ipaix)/ Jl^dx/ipa{x), with 
i^a{x) = [Da(x)]“. The second equality in Eq. d3 
follows immediately from x —x symmetry consider¬ 
ations. The ratios Nr/Nr and Tr/Tr have been com¬ 
puted numerically. The results plotted for a = 1/2 (Eig. 
El) and for a = 1 (Fig. E]) confirm the consistency of our 
phenomenological approach in both regimes. 

Clearly, our approach hinges on the phenomenological 
Eq. dSI and the explicit expressions we used for Va and 
Da{x). We now justify our choice for both. 


with P = P(x, (/, t), admits a uniform p.d.f., as one can 
prove by substitution; hence, the Po{x) of Eq. m with 
a = 1. The diffusion coefficient will be calculated again 
through Kubo’s formula: since in the stationary regime 
X and t are statistically independent, D = Di{x). Re¬ 
garding the drift velocity, the condition jo = 0 in Eq.Q 
requires that Va{x) = DQ,(a:)(ln[DQ,(x)Po(a:)])/ namely, 
for a = 1, vi(x) = D{(x), as expected. 

Coming back to the plots of Eigs. El and El we stress 
that: 

(i) The insets of Eigs.EKb) and (c) illustrate the depen¬ 
dence of Nr^r and Tr^r on the channel length L: Tr^r 
scale like whereas Nr^r grow insensitive to L. Of 
course, both statements are valid as long as L'^ 

_ (ii) Our expressions for TR,R{5a), adapted from Ref. 

[l3| . correctly reproduce the limiting values TR^R{6a) = 
L'^/SDq for Do <C Ds (gradient effects are superseded by 
thermal noise), and Tr,l(0) = L?/%Ds for Dq = 0 and 
Ja ^ 0 (purely active Brownian motion); 

(iii) On comparing the curves for Nr/Nr in Figs. EKb) 
andl^b) and those for Tr/Tr in Figs. E^c) andEKc), the 
different dependence of the two ratios on D^ at low ther¬ 
mal noise becomes apparent. This can be easily explained 
by inspecting the corresponding analytical expressions in 
the limit Do 0. For S^ = 0, > 0, i.e., 0 = 1/2, the 

two ratios are functions of Sy only and, therefore, inde¬ 
pendent of D^. For Sy = 0, S^ ^ 0, i.e., a = 1, they grow 
insensitive to uq, but do depend on S^ and, hence, D^j,. 
Accordingly, the limits Do ^ 0 and uq —>■ 00 coincide, as 
confirmed, for instance, by the numerical data in Figs. 
EKa) and (b). 
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FIG. 4: (Color online) Channel diffusion for = 0 and differ¬ 
ent 5^ and 77 ,^ (in legends): (a) Nr/Nl vs. vq for Do = 0.03; 
(b) Nr/Nl vs Do for vo = 1.5; and (c) Tl/Tr vs Do for 
Vo = 1.5. In (a)-(c) D^ = 0.1, L = 100/,^, and j/n = 5. 
The solid curves are the corresponding analytical predictions 
based on Eqs. ( 0 -® for a — 1. 


Finally, to fully answer our starting question (1), as 
the most likely exit end in the general case 7 ^ 0 and 
^ 0, we computed the rectification factor 




NlTr 


+ 1 


(9) 


where the ratios Nr/Nl and Tl/Tr are combined to¬ 
gether to quantify both the sign and magnitude of 
the symmetry-breaking mechanism responsible for the 
pseudo-chemotactic drifts. The most effective right (left) 
rectification corresponds to e = ± 1 , whereas for e = 0 
the opposite pulls by ^ 0 ( 2 ;) and D^{x) cancel each other. 


Note that for rjy = 1, = 2 and Sy = 5rj>, the transient 

rectification does not vanish [Fig. |2Ic)], even if under the 
same conditions Tr = Tr [Fig. [2{b), inset]. 


IV. CONCLUDING REMARKS 


The phenomenon of “drift without current” has been 
explained in Refs, [l^, [ 2 ^ using the phenomenological 
Eq. (|3]), by noticing that the statistical ensemble govern¬ 
ing the average drift (i.e., the rectification factor e in our 
notation) is different from the one required to compute 
the average current, ja = 0. The former consists of the 
representative points exiting an infinitesimally narrow 
neighborhood with coordinate x, with equal a:-dependent 
jump length in either direction, whereas the latter con¬ 
sists of all points crossing a channel cross-section with 
coordinate x at a given time, no matter what their jump¬ 
ing length. The two ensembles differ as an effect of mul¬ 
tiplicative noise [i.e., the x-dependence of Da{x)] and so 
do the currents thus calculated. 

Self-propelling artificial microswimmers reproduce 
that very same situation as a combined effect of nonequi¬ 
librium and the higher dimensionality of their dynamics, 
Eq. ( 0 ). In contrast to bacterial chemotaxis [l^, for 
an artificial microswimmer the self-propulsion parame¬ 
ters do not depend on the orientation. Here a depen¬ 
dence on the swimmer’s orientation might come into play 
due to, say, inertial or memory (i.e., nonlocal) effects, 
but surely not to some internal sensor-actuator path¬ 
ways, like in bacteria [I|. The microswimmers consid¬ 
ered here are characterized by very low Reynolds num¬ 
bers and small dimensions compared to the p{x) length 
scale; therefore, an orientation dependence of the swim¬ 
mer’s self-propulsion mechanism is not an option. 

For artificial microswimmers under the most common 
experimental conditions, is only weakly affected by 
the p gradient, while vo(x) is reported to grow linearly 
with p and then saturate at higher p P. [l4l - [l^ . The onset 
of “anti-Fick” cold-to-hot (pseudo-chemotactic) currents 
can thus be easily demonstrated. For instance, a source 
steadily releasing fuel into a JP suspension, causes a con¬ 
centration gradient around it; JP’s with x-independent 
rotational dynamics are driven away from the fuel source, 
whereas a tagged JP floating in such a depletion zone 
actually drifts toward the source. This prediction is in 
contrast with the experimental findings of Ref. Q , where 
Au-Pt micro-rods are reported to progressively cluster 
around an H 2 O 2 fuel source. If we assume that the self¬ 
propulsion model of Eqs. o holds good for a free swim¬ 
mer in the bulk (as established under the most diverse 
experimental conditions [^. [TTI. ) . the only explanation 
for such a discrepancy is that, upon migrating toward the 
fuel source, the JP’s come into contact with one another 
and eventually aggregate, as suggested, e.g., in Ref. [T7I[ . 
Drifts without current become observable at low swim¬ 
mer concentration. 
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